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We discuss the phase structure and the equation of state for QCD at non-zero temper- 
ature and density. Derivatives of In Z with respect to quark chemical potential fi q up to 
fourth order are calculated for 2-flavor QCD, enabling estimates of the pressure, quark num- 
ber density and associated susceptibilities as functions of fi q via a Taylor series expansion, 
^v^j . Also, the phase transition line for 2 and 3- flavor QCD and the critical endpoint in the (T, 

plane are investigated in the low density regime. 

\o ' 

. §1. Introduction 

(N 

We would like to report our recent study of QCD thermodynamics with a small 
but non-zero quark chemical potential \x q by numerical simulations. The interesting 
regime for heavy-ion collision experiments is rather low density regime, e.g. ^ q /T c ~ 
0.1 for RHIC and fJ- q /T c ~ 0.5 for SPS. The phase transition for 2-flavor QCD is 
known to be crossover at fj, q = and expected to become a first order phase transition 
at a critical endpoint. To find the endpoint is also an interesting topic in the low 
density regime, and it might be possible to detect the endpoint experimentally via 
event-by-event fluctuations in heavy- ion collisions. 

We discuss the equation of state at non-zero baryon number density in Sec. 2 ^ . 



The study of the equation of state gives the most basic information for the experi- 



ments. Quantitative calculations of thermodynamic quantities such as pressure and 
energy density are indispensable. In particular, since the number density fluctuation 
should be large around the critical endpoint, the susceptibility of the quark number 
density, given by the second derivative of pressure with respect to fj, q , is an important 
quantity. Several studies of the quark number susceptibility have been performed 
at fi q = 2 ). Moreover, measurements of pressure and energy density at \i q ^ 
were done using the reweighting method 3 ) , which allow the investigation of thermo- 
dynamic properties at non-zero baryon density. This approach, however, does not 
work for large fi q and large lattice size due to the sign problem. 

Our strategy is the following. We compute the derivatives of physical quantities 
with respect to /i g at \x q = 0, and determine the Taylor expansion coefficients in terms 
of /Uq 4 )' 5 ), in which the sign problem does not arise. Because the pressure is an even 
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Fig. 1. Coefficients of Taylor expansion, C2 (left) and C4 (right). To is T c at fi q — 0. 



function of fj, q , the //^-term is leading and /x^-term is the next to leading, and in fact 
only these two terms are non-zero in the high temperature (Stefan-Boltzmann) limit. 
We compute the Taylor expansion coefficients up to fourth order. The fourth order 
term enables us to evaluate the /^-dependence of the quark number susceptibility 
near \x q = . By estimating the change of the susceptibility, we discuss the possibility 
of the existence of the critical endpoint in the phase diagram of T and fi q ^ . 

In Sec. 3, we compute the Taylor expansion coefficients for the transition tem- 
perature T c 5 ) . Some investigations of the phase structure on the (T, fi) plane have 
been done by the rewighting method 6 ) and simulations with imaginary chemical po- 
tential 7 )' 8 ). Because the procedure to identify the phase transition point is complex, 
it is difficult to apply the method which we use in Sec. 2. Our method for T c is the 
following. We use the reweighting method, but we perform a Taylor expansion for 
the quark determinant and the fermionic operators and neglect unnecessary higher 
order terms of n q for the calculation of the derivatives which we want. Then, in 
the region where the sign problem is not serious, we calculate the Taylor expansion 
coefficients. This method reduces computer time drastically in comparison with the 
original reweighting method and enables us to use rather large lattice size. 

Finally, we discuss the critical endpoint again. Since a Taylor expansion must 
break down at the singular point, it seems to be difficult to find the position of the 
critical endpoint by the Taylor expansion. We focus on a property of 3-flavor QCD 
at very small quark mass 9 ). The crossover phase transition becomes a first order 
phase transition at a critical quark mass also for fj, q = 10 ) . Hence, by tracing out 
the position of the critical endpoint from that at ji q = to finite fj, q in the vicinity 
of jJLq = 0, we may be able to estimate the position of the critical endpoint for the 
real physical quark masses. We try it in Sec. 4. Section 5 presents our conclusions. 

§2. Taylor expansion of pressure in terms of n q 

Pressure is given in terms of the grand partition function Z(T, V, fj, q ) by p/T 4 = 
(l/VT 3 )\nZ. However, the direct calculation of \nZ is difficult, hence most of the 
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work done at fi q = for the calculation of pressure is done by using the integral 
method 11 )' 12 ^, where the first derivative of pressure is computed by simulations, and 
the pressure is obtained by integration along a suitable integral path. For fi q ^ 0, 
direct Monte Carlo simulation is not applicable; in this case we proceed by computing 
higher order derivatives of pressure with respect to /JL q /T at fj, g = 0, and then estimate 
p(Hq) using a Taylor expansion, 



p_ 

rp4 



T tflq 



+ £c n (T)(|n , (2d) 



Tfi n =l 



+ 6 



where c n = (l/n\)d n (p/T 4 )/d(fi q /T) n \ fMq=0 . Explicitly writing 

C2 = ^2, 04 = ^^(^4-3^), (2-2) 

„ / N f 9 2 (lndetM) \ . / / N { dflndetM) \ 2 

^2 - \— — gjp — ) + \yr — aji — ) 
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( Nf_ \ 3 <3 2 (ln dct M) ( 9(lndetM) \ 2 \ , / / Nf 9(lndctM) \ 4 \ 
1, 4 J 9^ V dft y 1 / \ V 4 9 M ) J ' 

for staggered type quarks on an iV 3 x A^ T lattice, and the coefficients for odd terms 
are zero. M is the quark matrix. \i is a chemical potential in lattice units, n = /i q a. 
Here, we used a property that the odd derivatives of In det M are purely imaginary 
and the even derivatives are real 5 ) . These derivatives are computed by the random 
noise method, which saves computer time. Furthermore, we do not need simulations 
at (T, n q ) = (0, 0) for the subtraction to normalize the value of p, since the derivatives 
of p\r=o,Li q =o with respect to fi q are, of course, zero. This also reduces computer time. 

We compute the pressure up to 0(/u, q ) using 2-flavors of p4- improved staggered 
fermion 13 ) at a bare quark mass ma = 0.1 on a 16 3 x 4 lattice. Then the quark 
number density n q and quark number susceptibility Xq are calculated up to 0(/U 3 ) 
and 0(fj, q ), respectively. These are obtained by the derivatives of pressure; n q /T 3 = 
d{p/T 4 )/d(n q /T), x q /T 2 = d 2 (p/T 4 )/d(n q /T) 2 , The details are given in Ref. 1). In 
Fig. 1, we plot the data for C2 (left) and C4 (right). Both of them are very small at low 
temperature and approach the Stefan-Boltzmann (SB) limit in the high temperature 
limit, as expected. The remarkable point is a strong peak of C4 around T c . 

We can understand this peak via two arguments. One is a prediction from the 
hadron resonance gas model 14 ) , which is an effective model of the free hadron gas in 
the low temperature phase. The model study predicts C4/C2 = 0.75 and our results 
are consistent with this prediction for T < T c ; in fact, as T increases 04/02 remains 
constant until T ~ T c , whereupon it approaches the SB limit. 

The other point is from a discussion of the convergence radius of the Taylor 
expansion. We expect that the crossover at ji q = changes to a first order transition 
at a point fi q /T c ~ 0(1) 6 )' 9 ). Then, the analysis by Taylor expansion must break 
down in that regime, i.e. the convergence radius should be smaller than the value 
of fi q /T at the critical endpoint. We define estimates for the convergence radius by 
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Fig. 2. Difference of pressure from \i q — (left) and Quark number susceptibility (right) as a 
function of T for each fixed fJ, q /T. Tq is T c at /i 9 = 0. 



Pn = \Z\ c n/c n +2\- We compute po and P2 from Co = p/T (0), C2 and C4. It is found 
that both po and p 2 are quite large at high temperature as expected from the SB 
limit, p'2 B — 2.01, pf B ~ 4.44. On the other hand, around T c , these are O(l), since 
C2 and C4 are of the same order, so that our results around T c suggest a singular 
point in the neighborhood of p q /T c = 1. 

Next, we calculate pressure and quark number susceptibility in a range of < 
Pq/T < 1 which is within the radius of convergence discussed above, Using the data 
of C2 and C4, p/T 4 and Xq/T 2 is obtained by A(p/T 4 ) = p(T, p q )/T A — p(T,0)/T 4 = 
c 2 (p q /T) 2 + c A (p q /T) 4 + 0(/^), and Xq /T 2 = 2c 2 + \2c A (p q /T) 2 + 0{p 4 ). 

We draw A(p/T 4 ) for each fixed p q /T in Fig. 2 (left) and find that the difference 
fromp| M(j= o is very small in the interesting regime for heavy-ion collisions, p q /T ~ 0.1 
(RHIC) and p q /T 0.5 (SPS), in comparison with the value at p q = 0, e.g. the 
SB value for 2-flavor QCD at p q = 0: p SB /T 4 ~ 4.06. The effect of non-zero quark 
density on pressure at p q /T = 0.1 is only 1%. Also, the result is qualitatively 
consistent with that of Ref. 3) obtained by the reweighting method. 

Figure 2 (right) is the result for Xq/T 2 at fixed p q /T. We find a pronounced 
peak for m q /T > 0.5, whereas Xq does not have a peak for p q = 0. This suggests 
the presence of a critical endpoint in the (T, p q ) plane. 

This discussion can be easily extended to the charge fluctuation xc- The spike of 
Xc & t T c is weaker than that of Xq-i which means the increase of the charge fluctuation 
is smaller than that of the number fluctuation as p q increases ^ . 



§3. Phase transition line in the (T, p q ) plane 



Next, we investigate the phase transition temperature as a function of \x q by 
performing a Taylor expansion around p q = O 5 ). In principle, the method via a 
Taylor expansion is applied for any physical observable. However, as the order of 
the derivative becomes higher, the number of terms increases exponentially and 
there exist serious subtractions, so the method in the previous section does not suit 
a quantity which needs complex procedures for the measurement such as T c . 
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We consider the following identity, 



(o) (M = (ow) (Aifi) /(w) iM , 



(3-1) 



where W = e W/4)(mdctM( M )-indctM(o)) e -S g (i3)+s g (f3o) for staggered type quarks, S g 



is the gauge action, and j3 = 6/g 2 . An expectation value (O) at (/3, fi) is computed 
by a simulation at (/3o,0). This is a basic formula of the reweighting method, but 
we expand lndetM(/x) and any operator of a fermionic observable using a Taylor 
series, e.g. lndet M{ji) - lndet M(0) = £~=i[<9(lndetM)/d/x](/i n /n!). This formula 
is valid in a region that the lndetM(^) does not have any singular point for each 
configuration. We moreover neglect the terms higher than 0(p, n ), then the resulting 
expectation value contains the error of the higher order of n but the error does 
not affect the calculation of the derivatives for the lower order than the neglected 
terms. In fact, if we expand into the form of the Taylor expansion from Eq.(3-1), wc 
obtain the correct coefficients of the Taylor expansion for (O) up to 0(/i n ), and this 
expression is much simpler than the method discussed in the previous section. Here, 
we try to trace out the phase transition line for non-zero /i using this method. 

However, if we want to use Eq.(3-1), we must discuss a famous problem, called 
"sign problem" . Because det M is complex at /x 7^ 0, if the complex phase fluctuates 
rapidly, the numerator and the denominator in RHS of Eq.(3-1) becomes vanishingly 
small, and then we cannot obtain the reliable answer for (O). Of course, the phase 
fluctuation is zero at /i = 0, but it becomes larger as /1 increases. However, in the 
context of the Taylor expansion, the applicable range of the reweighting method 
is rather easy to estimate by evaluating the complex phase fluctuation. For small 
/i = N^diq/T), the complex phase can be written by the odd terms of the Taylor 
expansion of In det M 5 ). Denoting detM = | detM|e ie , 



The first term is (A'f/4)Imtr[M~ 1 (<9M/(9/i)]/i. From these equations, we find explic- 
itly that the magnitude of 9 is proportional to /j, = N~ l (n q /T), V = and Nf. 
Moreover this can be computed by the noise method. Roughly speaking, the sign 
problem happens when the fluctuation of 6 becomes larger than 0{tt/2). Thus the 
sign problem is not serious for small \i and small volume, but that region becomes 
narrower and narrower as the volume increases. 

We observed the complex phase fluctuation on a 16 3 x 4 lattice 5 ) and found 
that the applicable range is not so small for this lattice size. The applicable range is 
narrower than the convergence radius discussed in Sec. 2, but it covers the interesting 
regime for the heavy-ion collisions at RHIC. We also found a tendency for the phase 
fluctuation to increase as quark mass decreases, and that the fluctuation is small in 
the high temperature phase. 

In the region where the sign problem is not serious, we determine the phase 
transition line. Because the first derivative is expected to be zero from the symmetry 
under exchange of to — fi, we calculate the second derivative of T c with respect 
to fi. Simulations are performed around the phase transition point (3 C by using the 




(3-2) 
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Fig. 3. Chiral susceptibility. Fig. 4. Sketch of the phase diagram. 



p4-improved action at ma = 0.1 and 0.2, corresponding m n /m p f» 0.70 and 0.85, 
on a 16 3 x 4 lattice. In Fig. 3, we plot chiral susceptibility as a function of j3 for 
H = 0, 0.05 and 0.1 at ma = 0.2. This calculation contains errors of 0(/U 4 ). From this 
figure, we find that the peak position of the susceptibility moves left as \i increases, 
which means that f3 c or T c decreases as fi increases. Assuming the peak position 
is (3 C , we determine the second derivative of (3 C . The truncation error of 0(/U 4 ) 
does not affect to this calculation. Combining with the results for the Polyakov 
loop susceptibility, we obtain d 2 f3 c /d/j 2 « — 1.1 and the quark mass dependence of 
d 2 /3 c /d/i 2 for ma = 0.1 and 0.2 is not visible within the accuracy of our calculation. 
(However, we find large mass dependence for very small quark mass. See Sec. 4.) 
The second derivative of T c is given by 




d 2 r c _ 

drf N 2 Tc dfJ 2 



(3-3) 



The beta function, a(df3/da), is obtained from the string tension data in Ref. 11). 
We then find (T c /2)(d 2 T c /d/i 2 ) = —0.07(3). We sketch the phase transition line from 
the curvature in Fig. 4. The diamond symbol is the critical point obtained by Ref. 6). 
The dotted line is the upper bound of the fit range to determine the curvature. At 
the relevant point for RHIC, this shift of T c is very small from that at ji = and the 
result is roughly consistent with those obtained by the other groups 6 )> 7 ). 

Comparing the result of the equation of state, we also find that the lines of 
constant pressure and energy density are roughly consistent with the phase transition 
line for small fj, q , which suggest the change of pressure and energy density is small 
along the transition line. Moreover, by measuring the Polyakov loop, we can observe 
that the external quark current is screened by dynamical anti-quarks 5 ^. 

§4. Critical endpoint for 3-flavor QCD 



Finally, we consider 3-flavor QCD. We perform simulations for Nf = 3 at ma = 
0.1 (m^/mp w 0.68) on 16 3 x 4 lattices, and at ma = 0.005 on 12 3 x 4 and 16 3 x 4 
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lattices, using the p4-improved action 9 ). The result of the curvature at ma = 0.1 is 
(T c /2)(d 2 T c /d^g) = —0.045(10), where we keep \i for strange quark equal to zero, 
and a(d(3/da) is obtained from the string tension data. The difference from Nf = 2 
is not large at ma = 0.1. However, a significant difference is observed between 
ma = 0.005 and 0.1 for Nf = 3. Using a perturbative beta-function, we obtain 
(T c /2)(d 2 T c /d / u2) = -0.025(6), and -0.114(46) for m = 0.1 and 0.005, respectively. 
Although the perturbative beta-function gives almost a factor of 2 smaller result 
than the non-perturbative analysis, this suggests that the curvature of the transition 
line becomes stronger as the quark mass decreases. 

The most interesting point for Nf = 3 is the existence of a critical quark mass 
m c on the fi q = axis, which separates a first order phase transition at small m 
and crossover at large m. In order to identify m c , we compute Binder cumulants 
constructed from the chiral condensate, B4 = ((Sip'ip) 4 } / ((Sipip) 2 ) 2 . It has been 
verified that the critical point belongs to the Ising universality class, i.e. the value 
of -B4 at T c for m c is the same with that of 3-dimensional Ising model 10 ) . Moreover, 
we expect such a critical point at fi q 7^ even in the region of large m, hence it is 
very interesting to investigate how the critical point moves as a function of jj, q . 

We apply the Taylor expanded reweighting method, explained in Sec. 3, with 
respect to quark mass, replacing fi by Am = m — mo, where mo is a simulation 
point, and we investigate the quark mass dependence around the simulation point, 
mo = 0.005. The reweighting factor is computed up to 0[(Am) 2 ]. This analysis 
suggests a critical value of m c = 0.0007(4) for fi q = 0. At this point, the Binder 
cumulant is consistent with the value of the Ising model, and also the peak height 
of the chiral susceptibility x^ shows scaling behavior indicating a first order phase 
transition, i.e. ~ V, f° r the volume size V = 12 3 and 16 3 below about this 

point. The corresponding value of pseudoscalar meson mass at the critical point 
is m£ rrf = 67(18) MeV. Hereafter, we use the value of m n at T = for setting a 
physical scale instead of the value of bare quark mass m. 

Next, we investigate the critical point for finite \i using the Taylor expanded 
reweighting method with respect to fj,, up to 0(fi 2 ). Since the sign problem is serious 
for the analysis on the 16 3 x 4 lattice, we analyze only the data on the 12 3 x 4 lattice. 
The data of B 4 suggests /x c "* = 0.074(13) or n° q rit /T = 0.296(52). An estimate for 
the transition temperature at corresponding is obtained from T c j sjo = 0.40(1) + 
0.039(4) {m^j^d) Xh \ Then the critical value is n° q rit = 52(10) MeV at the simulation 
point ma = 0.005. The corresponding m n is about 170 MeV. 

In the context of the Taylor expanded reweighting method for (2+l)-flavor QCD 
with equal up and down quark masses and a strange quark mass, m uc i and m s , the 
system is identical along the line of slope —2 in the (m U( i, m s ) plane near the m uc i = 
m s line. Along this line, the first terms of reweighting factor in terms of Am for m u d 
and m s are cancelled, i.e. if (m s — m^ 3 ^)/(m U( i — m ^ 3 ^) = — 2 is satisfied, the systems 
at i m ud, ms) and (m^ 3 ^,m^ 3 ^) are the same for m u d rs m s . We use this property for 
the extrapolation to the physical point. Using the lowest order chiral perturbation 
theory relation, m^/m 2 = (m u d + m s )/(2m U( i), this line in the (m u d,m s ) plane is 
translated to the line of slope —1/2 in the (m^m 2 ^) plane. Hence, when m 2 K = 
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—m\j1 + 3mi /2, the physics at (m^, m \) is roughly corresponding to that of 3- 

(3 ft 

flavor QCD with m\ . For the physical value m, = 140 MeV and txik ~ 500 MeV, 

(3 ft 

we obtain the corresponding pion mass for 3-flavor QCD, m\ = 416 MeV. The 
extrapolation of the critical surface to this pion mass value is performed by a straight 
line in the (m^,^) plane using the two critical points which we found at fi q = 
and /j,g ^ for 3-flavor QCD. We obtain an estimate for the critical point at the 
physical quark mass, ii c q nt k, 140 MeV. This value is smaller than the value of about 
240 MeV estimated in Ref. 6) for (2+l)-flavor QCD. Further study is clearly required 
to improve this analysis. 

§5. Conclusions 

Taylor expansion coefficients were calculated up to 0(/i ) for pressure and up 
to 0(/U 2 ) for phase transition temperature. We found that the differences of T c 
and pressure from those at fi q = are small for the interesting values of fi q for the 
heavy- ion collisions at RHIC and SPS. The fluctuations in the quark number density 
increase in the vicinity of T c and the susceptibilities start to develop a pronounced 
peak as \i q is increased. This suggests the presence of a critical endpoint in the (T, fj, q ) 
plane. Also, the critical endpoint for 3-flavor QCD near n q = was investigated. 
The extrapolation from the critical point at n q = suggests that the critical endpoint 
exists around fi q /T c ~ 1 for the physical quark mass. 
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